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BOUNDED ERROR SCHEMES FOR THE WAVE EQUATION ON COMPLEX DOMAINS* 

SAUL ABARBANELt, ADI DITKOWSKI*, AND AMIR YEFET* 

Abstract. This paper considers the application of the method of boundary penalty terms (“SAT”) to 
the numerical solution of the wave equation on complex shapes with Dirichlet boundary conditions. A theory 
is developed, in a semi- discrete setting, that allows the use of a Cartesian grid on complex geometries, yet 
maintains the order of accuracy with only a linear temporal error- bound. A numerical example, involving 
the solution of Maxwell’s equations inside a 2-D circular wave-guide demonstrates the efficacy of this method 
in comparison to others (e.g. the staggered Yee scheme) - we achieve a decrease of two orders of magnitude 
in the level of the ^ 2 -error. 

Key words. Maxwell’s equations, wave equation, finite difference time domain, error bounds, boundary 
conditions, complex geometries, staircasing 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. Hyperbolic systems of P.D.E.’s describing physical situations such as electro- magnetism, 
acoustics, elastic waves, etc, may under many circumstances be cast as wave equations for the various field 
components. 

One class of problems is that of solving numerically the Dirichlet problem on complex shapes, e.g., inside 
wave guides. For sufficiently non-simple geometries, the option of transforming the problem to body-fitted 
coordinates is not always a viable option, especially in three space dimensions. There are other options, such 
as using Cartesian grids and approximating the body shape via “staircasing” , “diagonal split cell model”, 
etc (see for example Chapter 10 in reference [4]). It is well known that these devices are not very efficacious, 
particularly in the high frequency regime. We shall demonstrate that “staircasing” can fail even for low 
frequencies. 

In this paper we consider the application of the method of boundary penalty terms (“SAT” , see references 
[1], [2], [3]) to the numerical solution of the wave equation in a finite domain with Dirichlet boundary 
conditions. 

In Section 2 we develop the theory that allows us to use a Cartesian grid on complex geometries and yet 
maintain the order accuracy with a linear temporal error-bound. 

In Section 3 we construct a second order accurate scheme that fulfills the conditions imposed by the 
theory presented in Section 2. 

Section 4 is devoted to a numerical example the solution of the transverse magnetic (TM) Maxwell’s 
equations [4] between two concentric circles. (This configuration might be considered as a cross-section of a 
very long wave-guide.) This problem is solved using four different numerical algorithms. Two of them solve 
the first order system with “staircasing” - the Yee staggered scheme [6] and a 4 th order spatially staggered 
scheme due to Turkel and Yefet [5]. The other two solve the wave equation directly on a non-staggered 
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Cartesian grid, one with the SAT formulation and one without. All three “standard” (non-SAT) algorithms 
have very large errors; the SAT algorithm has errors that ar? at least two order of magnitude smaller. 
Summary and conclusions, and ideas for future work are presented in Section 5. 

2. Theoretical Framework of the Method. In reference [1], [2] and [3], it was shown how the case 
of a one-dimensional P.D.E. can be used as a building block for the multidimensional case for constructing 
error-bounded algorithms over complex geometries with Dirichlet boundary condition. We therefore start 
with the following one- dimensional problem: 

d 2 u d^u 

(2.1) + f(x, t); T l <x<T r , t> 0 

(2.0a) u(x, 0) = uo(x) 

(2.0b) — u(x,0) = u to (x) 

(2.0c) u(F L ,t) = g L (t) 

(2-0d) u(T R ,t) = g R (t) 

and f(x, t) £ C 2 . 

Let us discretize (2.1) spatially on the uniform grid presented in Figure 2.1. Note that the boundary 
points do not necessarily coincide with X\ and x at- Set Xj + 1 - Xj = h, 1 < j < N — 1; x\ — = 7 

0 < 1l < 1; T r- xn — liih, 0 < 7h < 1. 



Fig. 2.1. One dimensional grid. 


Since, unlike the cases discussed in [1], [2], equation (2.1) has a second time derivative, attempts to 
apply naively the methods presented there fail. The reason is tkat if we follow the procedure used there and 
write the following discrete approximation to (2.1), 

(2.2) ^u = Du + f(i) + ? e 

where u is the projection of the exact solution u(x, t) onto the grid, i.c. u(xj, t) = Uj(t) = u(£); and write 
the numerical scheme 

(2.3) ^7 = [D\ - t l (A l v - g L ) - t r (A f v - g ft )] + f (t) , 
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then the equation for the error vector e = u — v becomes 

(2.4) + 

In the above, v is the numerical approximation to u, and 

(2.5) M = D — — trAr . 

D is a differentiation matrix of the proper order of accuracy that does not use boundary values. The matrices 
Ai and Ar arc defined by the relations 

(2-6) Alu = gL ~ T l, = g/? - T r , 

i.e., each row in A^(Ar) is composed of the coefficients extrapolating u to its boundary value gL{gii) at 
^l{^r) to within the order of accuracy. (The error is then Tl(Tr).) The diagonal matrices tl and r R are 
given by 


t l = diag (TL 1 ,n, 2 ,-*-,TL JV ); r R = diagfr^ , tr 2 , • • • , t Rn ) . 

The constrain on the construction of the A 1 s, r’ s and D is that M in (2.4) be negative definite. The negative 
definiteness of M is a necessary condition for extending the 1-D theory to the multidimensional case (see 
[1] , [3] ) . Also in (2.4) 


(2-7) T = T e - t l T l - t r T r = (7\, T 2 , ■ ■ ■ , T m , ■ ■ • , T n ) t . 

If the matrix M can be diagonalized* , then 

(2.8) M = Q~ l KQ 


with the diagonal matrix, A, having the eigenvalues of M. Defining /i = Qc, equation (2.4) becomes 


(2.9) 


d 2 fx 
dt 2 


= A fx + QT 


= A /i + T . 


This is an un-coupled system of O.D.E’s. The general solution for the m th equation is: 

^ rt 

A^m(^) = ^ m T c m2 e ^ 771 + .. I sinh ( \J A m (t s)^T 1 m (^)ds . 

\ *rn 

Recalling that at t = 0, e = e* = 0 (i.e. ja — fx t = 0 at t = 0), the solution of (2.9) becomes: 

i 

(2.10) Hm(t) = -7= / T m {s) sinh [\/A^(l - «)]ds ■ 

\Mm 7o 

Note that unless all the eigenvalues of M are real and non-positive some of the ^/A m ’s will have a positive 
real part, in which that case at least one of the /i m ! s may grow exponentially in time. In order to prevent 
this, we have to demand that M, in addition to being negative definite, also possess only real eigenvalues. 

•Extensive numerical evidence has shown that the M in [1],[2] (i.e. representing the second derivative to 4 th and 2 nd order 
accuracy, respectively) has distinct eigenvalues and hence is diagonalizable. 
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Furthermore, in order to use the 1-D scheme as a building bio ok for multidimensional schemes, M should 
be built in a way that verifies that the property of real negative eigenvalues carries over to the multi- 
dimensional differentiating matrix. One way to achieve this goal is to construct M as a negative-definite 
symmetric matrix. Then an estimate on the error bound can be derived directly from the solution (2.10), 


|Mm(0l — 




-T t 

7 TTlAr L 


where T mM — maxo< s <( |T m (s)|. Then, for a normalized Q, 

(2-11) ||6|| = ||M||<l||f M p, 

Co 

where Co = minm-i,...^ Therefore || e\\ grows at most linearly with t. 

This result, of a linear temporal bound on the error-norm, can also be derived by resorting to energy 
method (see [3]), instead of directly from the solution. 

Also, as mentioned before, the construction of multi- dimensional case 


d 2 u 
dt 2 


V 2 u + /(x, t] 


on complex shapes is completely analogous to the method indicated in [1], [3]. 

3. Construction of the Scheme. This section is devoted to the task of constructing a symmetric 
negative definite matrix M for the case of a second order accurate finite difference algorithm. 

Let 


r 

~ 1 -2 1 0 


1-210 


0 1-21 


0 0 1-21 


1-2 1 00 


1-2 10 


1 -2 1 


1 -2 1 


’ 0 


i — 

o 

o 

o 

0 

1 

C2 


1-3 3-1 

C3 


-14-6 4 -1 

CN-2 


-14-6 4 -1 

cat-1 


-1 3-3 1 

0 


0 0 0 0 
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(3.1) 


where 

(3.2) 

and 


I 

0 

0 

0 

1 

> 

0 1-21 


-12 0-21 


-1 2 0-21 


-1 2-10 


0 0 0 0 

j 


_ . cat-i - c 2 

Cfc — c 2 -\ — — g ( k - 2) , 


(3.3) 


c = 


C/y^i — f>2 

A T -3 


Note, that as in [2] and [3], we had to resort to using connectivity terms , the last two matrices in (3.1). 


(3.4) 


A, = 


j( 2 + 7l)(l + 7 l) ~7l(2 + 7l) 2(7l + 7l) 0 ••• 0 


l(2+7£.)(l + 7i) — 7 z.(2 + 7l) 1(7l+7l) 0 ••• 0 


(3.5) 




0 ... 0 -(7fl+7fl) -7fl( 2 +7«) + + 7«) 


0 ••• 0 + -7/?(2 + 7r) -(2 + 7ii)(l + 7n) 


(3.6) 


n = Jo diag [tl, , tl 2 , rt 3 , 0 , . . . , 0, 0] ; 


(3-7) t r = 7^ diag [0, 0 , . . . , 0 , t R n _ 2 , t R n _ 1 , r^] ; 

In order to make the matrix M = D — tlAl — trAr symmetric we choose: 


c 2 = 


(1 -'ll,) 7 l 
2 


Cn~ 1 


7~[ j 2 


(1 ~ 7 r) 7 R 
2 

3-7 l - 27 L r Ll 

1 + 7L 
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(3.8) 


(3.9) 




TRn-1 


T Rn- 2 


-2 + 7 l + 7 L tli 

2 + 7l 

3 — 7k — 2 7r r Rw 

1 + 7fi 

-2 + 7w + 7/» Tfiiv 

2 + 7fl 


tl 1 ,t Rn > 4 . 


The proof that the symmetric matrix M is indeed negative-definite is given in the Appendix to this 
paper. 

Note also that instead of solving (2.3) directly as a 2 nd order O.D.E. system in time, one can solve 

c/w 

— = [Dv - t l {A l v - gi.) - t r (A r \ - gfi)] + f 
at 

dv 

— = w . 
dt 

(3.10) 


The number of ’variables’ has increased from N to 2 TV but one gains in the simplicity of the time integration. 

4. Numerical Example. Wc consider the dimensionless Maxwell’s equation for transverse magnetic 
field (TM, see [4]) in two space dimensions: 


(4.1) 

(4.2) 

(4.3) 


d£ _ dHy _ dth 
dt dx dy 
dH x _ dE 
dt dy 

dH y _ dE 
dt dx 


where H x and H y are the x and y components of the magnetic vector, H, and E is the electric field in the 
z-dircction. The set (4.1) (4,3) is to be solved in the space between two concentric circles, g < r < 5. We 
consider the case of perfectly conducting boundaries. Thus the Doundary conditions are given by 


(4.4) E{\,e,t) = 0 

(4.5) E(l,0,t) = O. 

We choose the following initial conditions (note the polar coordinates r,0): 


(4.6) 


(4.7) 


(4.8) 


E(r, 0 , 0) = cos 6 [J\ (uor) + a Y\ (u;r)] 

H y (r , 0, 0) = — sin 2 o{ [ J\ (u?r) + a Y (u;r)] 

L 2 ujt 

[Jo(c^r) - J2(ur) + a Y 0 ur) - a T 2 (^)] j 
cos^ 0 

H x (r, 0, 0) = [ Ji (wr ) 4- a Yi («r)] 

ujr 

. 2 

— [J 0 (^r) - Ji (wr) + 1 y 0 (wr) - a ^(wr)] 


where the J n ’s and the 


T n ’s are Bessel functions of the first anc- second kind of order n, respectively. Also, 


(4.9) 


a 9* 1.76368380110927; u ^ 9.813695999428405 . 
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The exact solution of the IBV problem (4.1)-(4.8) is given by: 


(4.10) 

(4.11) 

(4.12) 


E(r,0,t) = cos (cot + 8) [ J\ (u;r) + a Yi (cor)} 

H y (r y 8,t) = — — cos 0 cos(c<^ + 0) [Ji (cor) -f a Yi(cor)] 
cor 

-h| cos0sin(a;t + 0) [Jo(wr) — J 2 (cor) -f a Y$(cor) — a Y 2 (cor)} 
H x (r, 0, t) — — - cos 0 cos (cot + 8)[J\ (cor) + a Y\ (u;r)] 

— \ sin0sin(u;t + 8) [ Jo(cor ) — J 2 (cor) + a >o(wr) — a Yz(cor)} 


We note that we can extract from (4.1)- (4.3) a wave equation for the electric field E, 

/A d 2 E d 2 E d 2 E 

( ' ) dt 2 dx 2 + dy 2 ' 

The boundary conditions on E in (4.13) arc given by (4.4) (4.5). The initial condition E(r , 0 , 0) is given by 
(4.6). We need an additional initial condition on E t , which we obtain by differentiating (4.10), namely 

(4.14) Et(r,0, 0) = — co sin 0 [ J\ (cor) + a Yi(cor)\ . 


Four numerical schemes were used to solve the problem: 

(i) The Yee scheme [6]. This second order accurate scheme is staggered both in space and time. This 
entails putting initial conditions of H x and H y at At/2 rather than at t = 0. These initial conditions 
are derived from the exact solution. The numerical solution is carried out on the “staircased 11 domain 
shown in Figure 4.1. 

(ii) A modification of the Yee scheme (designated Ty(2 ; 4)), sec [5]. This one has 4 th order spatial 
accuracy and 2 nd order in time. The stagger and the “staircased” domain are maintained as before. 

(iii) The SAT algorithm for the wave equation described in Sections 2 and 3. The grid used for the 
numerical integration is shown in the right side of Figure 4.1. The time evolution is done by a 4 th 
order Runge-Kutta method. 

(iv) An algorithm which formally looks like the SAT in (iii), but is applied to the “staircased” domain 
of Figure 4.1 (rather than SAT one). To order 0(h 2 ), this is equivalent to using a standard spatial 
central differencing scheme with the nodal points at edges of the domain given the boundary value 
zero. The time integration is done as in the case (iii). 

We first present the L 2 error in E for all four schemes at t — 1 and t = 10 for the cases Ax = Ay = h = 
1/40, h = 1/80 and h = 1/160, see Table 1. At was 2/3 h for the Yee scheme, h/ 18 for the Ty(2,4) scheme 
and h/ 5 for the SAT schemes. 

It is immediately apparent from the table that the SAT-error (scheme iii) is at least 2 orders of magnitude 
smaller than that of the other three algorithms at all the various times and grid spacings. 

Since the non- SAT schemes have errors which are unacceptably large we do not show details of their 
temporal behavior. The SAT algorithm (scheme iii) has an L 2 error which grows in time as shown in Figure 
4.2. We see that this temporal growth is bound by a linear curve, whose slope depends on h. We note that 
for all reasonable dimensionless time the error is quite small, especially for h < 1/80. 

5. Conclusions and Discussion. 

(i) It seems quite clear from the evidence that the failure of the non-SAT schemes is due to the fact 
that “staircasing” misrepresents the shape of the body. In the SAT scheme, on the other hand, the 
penalty terms take account of the true shape. 
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Fig. 4.1. The “staircased” domain (left) and the SAT grid (right), h = 1/40. 




h= 1/40 

h= 1/80 

1/160 


t = 1 




i 

Yee 

0.4322 

0.3635 

0.1742 

ii 

Ty(2,4) 

0.4038 

0.3347 

0.1579 

iii 

SAT 

0.001203 

0.0001705 

1.5019e-05 

iv 

Staircased 

0.1022 

0.05041 

0.01936 



h = 1/40 

h = 1/50 

h= 1/160 


t = 10 




i 

Yee 

0.5101 

0.4364 

0.6683 

ii 

Ty(2,4) 

0.2642 

0.7079 

0.7243 

iii 

SAT 

0.008435 

0.0008354 

8.2707e-05 

iv 

Staircased 

0.7929 

0.473-5 

0.7829 


Table 4.1 
The L% error. 


(ii) The numerical results validate the theoretical predictions of the temporal behavior of the L 2 norm 
of the error. 

(iii) Grosso-modo the CPU time per node is of the same ore er for all schemes. 

(iv) The results from Table 1 and Figure 4.2 seem to indi:ate that the scheme (iii) converges as /i 3 , 
although the algorithm has a truncation error of order /i 2 . We do not understand this pleasant 
anomaly, although it is possible that even with h = 1 / 1 30 we are not yet in the asymptotic conver- 
gence regime. 

(v) In the future, we would like to apply the SAT methodology directly to hyperbolic systems such as 
(4.1) (4.3). The theory is not complete yet. 
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L2 error 



Fig. 4.2. SAT, L2 error vs. time. 


Appendix. We decompose the matrix M, defined in (2.5) and (3.1) to (3.8) as follows: 


(5.1) M — — [otM 1 + (1 — a)M 2 ■+■ 4- A/4 + A/5] 

where: 


(5.2) 


-2 1 
1 —2 1 

1 -2 1 


Mi = 


1 


-2 1 
1 -2 1 

1 -2 


(5.3) 


M 2 


0 0 0 

0 0 0 

0 0-1 1 
1 -2 1 

1 -2 1 
0 1-100 
0 0 0 

0 0 0 
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r 

m 4 

1,2 

m 4 

1,3 

m 4 




1,2 

m 4 

2,2 

m 4 

2,3 

m 4 

0 

( 5 . 5 ) 

m 4 = 

1,3 

m 4 

2,3 

rn 4 

3,3 

m 4 ’ 





0 


0 


where: 
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0 


o 


(5.6) 

M s = 

iV— 2, JV — 2 
m 5 

N-l,N~2 

m 5 

N,N — 2 
m 5 



H _N-l,AT-2 
U ra 5 

/V-l./V-l 

m 5 

Af,N-l 

m 5 



N,N—2 

m 5 

N,N- 1 
m 5 

N,N 

m 5 


where: 


m"' N = l + 2a- (1 + 7fl) + 

5 2 

m w,JV-i = _ 2 - a + 7 h (2 + 7h) t Rn 

m N,N-2 _ , _ 7* (1 + 7ft) 

5 2 


/V-l.AT-1 


= 2a + 


7 7fl ~ 4 (l + t|) ~ 7ft 2 ( 2 + 7ft) C 1 + 4 tr„) 

2(1 + 7*) 


JV-l.AT-2 , 37 r 7r 2 2 


Tfl 


N — 2,N —2 


= 2 a + 


-4 + 7ft 2 - 7 r 3 - 7r 2 (1 + 7r) tr a 

2 (2 + 7r ) 


The matrix M\ is negative-definite and bounded away from 0 by h 2 7 r 2 by the argument leading to eq. 
(2.4.31), see appendix to chapter 2 in [3]. A /2 is non-positive definite, sec eq. (2.4.34) and (2.4.35) in that 
appendix. From (3.2), (3.3) and (3.8) follows that > 0, k = 1 , . . . , TV, therefore, the matrix A /3 is non- 
positive. For a given value of 0 < a < 1 , tl, 1 and tr n can be found such that the matrices A /4 and A /5 will 
be non-positive, for all and y R . For example: for a = 1/10, t Li = = 4; for a = 1/2, r Ll = = 9 

and for a = 8 / 10 , = 24. This completes the proof that M is indeed a negative- definite matrix, 

bounded away from 0 by a 7 r 2 . Therefore the norm of the error vector || e || can grow at most linearly in 
time, see equation ( 2 . 11 ). 
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